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Abstract 

Modified gravity provides a possible explanation for the currently observed cosmic accelaration. 
In this paper, we study general classes of modified gravity models. The Einstein-Hilbert action is 
modified by using general functions of the Ricci and the Gauss-Bonnet scalars, both in the metric 
and in the Palatini formalisms. We do not use an explicit form for the functions, but a general 
form with a valid Taylor expansion up to second order about redshift zero in the Riemann-scalars. 
The coefficients of this expansion are then reconstructed via the cosmic expansion history mea- 
sured using current cosmological observations. These are the quantities of interest for theoretical 
considerations relating to ghosts and instabilities. We find that current data provide interesting 
constraints on the coefficients. The next-generation dark energy surveys should shrink the allowed 
parameter space for modifed gravity models quite dramatically. 



I. INTRODUCTION 



The discovery of the current accelerated expansion of the Universe has generated a 
lot of excitement in the last few years. It turns out that this unexpected behaviour can be 
modeled in many ways, most simply through a cosmo logical constant. On the other hand it 
has been argued that, since this behaviour appears on large scales, the acceleration of the 
universe is due to some modification of gravity on such scales. 

There are quite a few models which can describe such modifications of gravity, such 
as scalar-tensor theories Q]-[l5], the so called f(R) theories 16]- 
scalar tensor theories), brane world models of which Randall-Sundrum models 
models 



26j (which are a subset of 
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30 



271 and DGP 



311 ] are specia 
invariants, such as f(R,R,Q B ) [32 



cases, and some more complicated function of curvature 
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36j where Rq B represents the Gauss-Bonnet 



combination defined as R 2 - AR aP R af3 + R a ^ 5 R a ^ s . 

Modifying gravity in a consistent way is not an easy task, hence many of these models face 
quite stringent theoretical bounds which can reduce the parameter space of these theories. 
In the same way one can try to use data to restrict further the parameter space and possibly 
rule out some of these models. 

In this paper we study f(R) and f(R,RQ B ) models, most generally, without specifying 
any explicit form for these functions. In the f{R) case, we consider both the metric based, 
and the Palatini formalisms. We use the recent expansion history of the universe, recon- 
structed allowing the Hubble parameter to be a free function in redshift bins, using type la 
supernovae (SNe la) data together with relevant contraints from the cosmic microwave back- 
ground (CMB) and galaxy surveys, and invert this to obtain constraints on the parameters 
of the f(R) models. This perspective has not been adopted before. Although it is similar 
in nature to the reconstruction of the quintessence potential 3?j], here the scalar degrees of 
freedom come entirely from the gravity sector. 

Much work has been dedicated to the subject of constraining f(R) theories using solar 
system measurements 38(- 44J. We will not study this here for several reasons. First, as be- 
lieved by some (e.g. 45|), the cosmological parameters may be different from those measured 
locally, because as the background changes from a locally spherically symmetric metric to 
a homogeneous and isotropic one, the behaviour of quantities such as the coefficients of a 
Taylor-expansion of f{R) might be quite different. Second, for this theory locally (at the 
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solar system scale) the weak field approximation may not hold, so that it cannot be matched 
to a perturbation of GR; i.e., the usual constraints cannot be trivially applied because these 
are found by assuming the metric to be a perturbation about the GR-Schwarzschild solution 
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Furthermore we are restricting our study to the background evolution at low redshifts, 
assuming GR-like behaviour until then. We are not taking into account perturbations and 
their evolution, nor considering the stability of modes in the early universe, in constraining 
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the models. These relevant issues have been studied in 

This paper is organized as follows. In Section II we discuss the general metric based f(R) 
case and its solutions. In Section III we discuss the Palatini formalism and its solutions. In 
Section IV we consider the generalization including the Gauss-Bonnet term. We end with a 
discussion and conclusions section. Appendices follow. 



II. GENERAL METRIC-BASED f(R) MODEL 

Let us begin with an f(R) theory with the following action 

S= f^^^ + S m , (1) 
J lb7rG 

where R is the Ricci scalar and S m is the action for the matter fields. The equations of 
motion are 

(1 + fn)G, u - l -g^U - Rf R ) + g^Df R - V,VJ R = SttGT^, 

where an underscore R implies a partial derivative with respect to R (f R = df /OR). In a 
FRW background then, 

3(1 + f R )H 2 + \{f~ Rf R ) + 3H 2 f' R = 8nGp , (2) 

where H is the Hubble parameter and a prime denotes differentiation with respect to N = 
ln(a/a ). Evaluated today this equation becomes 

l + /? + i[«-6/3(fl 7 /ffo + 2)]+ 7 ^ = ^Pd, (3) 

where an underscore implies present values, and a = fo/H 2 , (3 = f R0 , and 7 = f R mHl 
are the dimension-less parameters of this theory. Since we have three parameters to solve 
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for, we need three equations. These are obtained by differentiating the Friedmann equation 
twice. The equations then involve higher derivatives of H evaluated today (in this case upto 
four) . 

Assuming that / can be Taylor-expanded in R about today and retaining only up to the 
second order terms (this being the simplest non-trivial case), 



/ « fo + fm Hi 



R — Rq 



i 

2 



+ fRRo H 



R — Rq 



(4) 



In order to use the equations of motion effectively one needs to truncate the Taylor expan- 
sion at some order, otherwise one would need an infinite number of equations, obtained from 
repeatedly differentiating the Friedmann equation, to specify an infinite number of param- 
eters, the Taylor coefficients. We chose to truncate this expansion of f(R) at second order 
because stability constraints on these theories involve f R and f RR terms (see e.g. j^). 
Differentiating the Friedmann equation once gives 



6 HH' (1 + f R ) + 3 (HH> - H 2 ) f' R + 3 H 2 f' R = 8nG p' 



(5) 



This relation together with the conservation of stress-energy (which ignoring contribution 
from radiation leads to) 

p' = -3(1 + -3p . (6) 

is equivalent to the second Einstein equation. Equations fflh, (jHJ) and (jSJ) lead to 



6 W (1 + fm) + 3 {W - 1) f RR0 R' + 3 fjua R' ' = -9 Q r , 



(7) 



where we have defined the (present day) derivatives of the (normalized) Hubble parameter 

as 

H' = ^ and H" = ^f. (8) 

Hq Uq 

For a flat Friedmann-Robertson- Walker background 



R = 6(HH' + 2H 2 ) 

R' = 6 (H' 2 + H H" + AH H') 

R" = 6 (H H'" + 3H' H" + 4 H' 2 + AH H") , 

hence equation (j7j) can also be written as 

2W + 3fi m + 6 7 (H m + 3H" + AH' H" + H' 3 + 7H' 2 - AH') 



(9) 







(10) 



Differentiating the Friedmann equation a second time, evaluating it today and using 
equation ffTUl) gives 
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\ fl m (3W + H' 2 + H") x [m' A + 15H' 2 H" + 6H' 3 H" 

- m" 2 + m"H"' + 6W V" - n"n"' + w'w" 



where 



H" 



TJIII 



and 



H" 



HIT 
Ho 



We see that 



a = a(H',H",H'",H m ,tt r 



(11) 



(12) 



(13) 



and the same is true for (3 and 7. The relations are non-linear. 

Posterior distributions of TC', TC", Tl'", TL"" and fl m are obtained from a likelihood analysis 
of data in the following way. We use 182 Type la supernovae (SNe la) from the HST/GOODS 



program, together with first year SNLS and some nearby SNe la, as 



(R, l a , Q b h 2 ) combination, where R and L are CMB shift parameters [5 ll]. is used to account 



for relevant constraints from the CMB 



compiled by 



50|. The 



521 ] . The SDSS baryon acoustic oscillation (BAO) 



scale measurement is also used 53|]. Following the analysis method of 54J- 56j], 5jj, model 
independent constraints are derived on the Hubble parameter in linear redshift bins using a 
Markov Chain Monte Carlo (MCMC) algorithm. The if (z)'s of the MCMC chain elements 
are then converted into the derivatives of 7i, as described in Appendix A. The top panel of 
Fig 1 shows the constraints thus derived on what we for convenience shall sometimes refer 
to as the "initial parameters". From the MCMC chains now we can obtain constraints on 
a, P, and 7 using the equations derived above. 



Linearization 

Equation (1111) is highly non-linear in the initial parameters. Small changes in these 
parameters can have a large effect on the values of the f(R) parameters, so that a relatively 
small number of high likelihood combinations of the initial parameters can lead to a large 
range of values for the f(R) parameters. This leads to flat ID parameter likelihood curves. 
This effect is suppressed in the MCMC posterior of the f(R) parameters, because in addition 
to the likelihood this takes into account the number of samples that fall into each bin. In 
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FIG. 1: Results for general metric-based f(R) gravity: constraints on a = fo/H^, [3 = /ro, and 
7 = IrroHq (bottom panel), obtained from constraints on £l m and the derivatives of H (top 
panel), using current cosmological data. 

other words the discrepancy between these curves is due to there being a main posterior 
peak (representing a concentration of high likelihood points) together with isolated points in 
parameter space that are allowed by the likelihood. In addition, some allowed values of the 
initial parameters lead to singularities; these are the zero values of the 7i"s and the zeroes 
of denominator in the expression for 7. Given the discrete samples in the MCMC chain, the 
parameters a, (5 and 7 don't actually blow up, but instead as a result have a larger range. 
These problems should be at least somewhat eased by better data. For now, in order to 
avoid the discrepancy between the likelihood and the posterior, and to discount the region 
around singular points, we proceed to find solutions by first linearizing the equations for the 
modified gravity parameters about the mean values of the initial parameters, supported also 
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FIG. 2: 13 as a function of H'" and ft"" setting ft' and ft" to their best values. 

by the fact that the initial parameters have close to Gaussian distributions. This approach 
will be used in subsequent sections as well. 

We linearize the equations for a, f3 and 7 about their maximum likelihood (in MCMC, 
the relevant quantity is the mean) values, which are obtained in turn from the mean values 
of the initial parameters: 



Pj=Pj 



Pj=Pj 



(14) 
(15) 
(16) 



da 

a « a + {pi - pi) — 
opt 

(3 w $ + ( Pl - p.) |^ 

OPi 

7 ~ 7 + {Pt ~ Pi) -K- 

°Pi p j=Pj 

where pi = (TL f , Tt", TC'", TC"", fi m ) and hats represent mean values. 
We obtain 

a = -29.7708 + 31.9923 H' - 0.652908 H" - 0.0628292 H'" 

+ 0.0513295 H"" + 63.1444 fi m , (17) 

(3 = -1.89148 + 2. 29135 ft' - 0.0479787ft" - 0.00179758 ft"' 

+ 0.00402273 ft"" + 4.20931 Vt m , (18) 

7 = -0.00179067- 0.00887807ft' + 0.000625996 ft" + 0.0000797733 ft'" 

- 0.0000138531 ft"" - 0.02083 VL m . (19) 



Even though the linearization is about the mean values of the initial parameters, we 
then use the MCMC chain to obtain the distributions of the modified gravity parameters. 
These are shown in the bottom panel of Fig 1. As described above these represent the 
main body of solutions of the f(R) theory allowed by the data; given the non-linearity of 
the equations there are other solutions isolated in parameter space. We see that while the 
order of magnitude of a and ft is constrained, 7 oc fm^o is found to be slightly negative 
today over most of its allowed range. Even though we are not discussing stability issues 
here, 
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24 



26| find that Jrr at high redshift needs to be positive in order to avoid 



instability and obey GR. Therefore under this theory 7 would have to change sign at some 
R before today. This is interesting. 

Fig 2 shows the solutions for ft at the mean values of the better constrained derivatives 
TC' and Tl", and over the allowed ranges of the higher derivatives. A range of values for ft 
are possible, including the special case of ft — 1/3 discussed below. 



A. Metric based f(R) with ft = 1/3 

If we assume that in the solar system the real metric can be expanded about GR 

Schwarzchild, then for the previous action, the effective Newton's constant can be writ- 
□ 

ten as [12J 

G 4 ° 



3 1 + /*' 

which can be recast as a constraint on ft today (imposing G = Geff,o) 

fm = ft=\- (20) 
Imposing ft = 1/3, equation d3J) and (JTU1) give 

\ a - \{H' + 2) + 6 7 (H" + H' 2 + AH') = Q m - f , (21) 

and 

§ W + 67 3H" + W (h' 2 + m' + AH" - 4) + W" =-3fi m . (22) 

The data define constraints on the H n s as discussed in the previous section. Equations ([2"T]) 
and (1221) can then be used to solve for a and 7. The complications arising due to the 
equations being non-linear in the initial parameters, as discussed earlier, apply here as well. 
Hence as before we proceed to linearize the equations of motion about the mean values of 
the TC n s and Q m . 

8 



Linearization 

Using equations (|T4|) and (|T6|) together with equations (I2TI) and (I22l) . we get 

a w -4.53219+ 12.1526 W- 1.54509 W- 0.974158 W" + 20.5967fi m , (23) 
7 w -0.028308 - 0.0521268 ft' + 0.0262386 ft" + 0.00930701 ft'" - 0.1394561^ . (24) 

These relations, used on the MCMC chains, give the results shown in Fig 3. 
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FIG. 3: Same as Fig 1 for the j3 = g case of the general metric-based f(R) theories. 

III. f(R) IN THE PALATINI FORMALISM 

So far we have considered theories for which the gravity variables were chosen to be 
the metric elements g^ u . On the other hand, after writing down the action, one is free to 
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choose a different set of fields. For example, one can choose to find the equations of motion 
by varying the action with respect to the following two tensorial quantities, 5g^ v and ST^ U , 
which are the metric and the Christoffel symbols perturbations respectively. In standard GR, 
the approach of choosing two different field variables, leads to the same standard Einstein 
equations of motion. However, in f(R) theories, with f R ^ constant, the equations of motion 
are indeed different. Of course, one may wonder which description of gravity is the correct 
one. However, since gravity is the least well known force, many physicists have argued that 
the possibility of introducing extra fields should be considered. 



In this formalism (refered to as Palatini, also studied in [19|, |57|, |58|, [59] ) , R pu becomes a 
function only of the Christoffel symbols, so that it cannot be written as usual in terms of 
derivatives of the scale factor. In this case it is customary to write the general action in the 
following way 

S= [d 4 xV^g{^ + S m . (25) 
J 1o7tG 

The Einstein equations obtained by varying the action with respect to the metric elements 
are 

f R R^-\fg^ = ^GT iaf , (26) 

the trace of which is 

f R R-2f = -8nGp m . (27) 
Taylor-expanding / in terms of R about today, as before, we have 

f=f=a + P(R-Ro) + kl(R-R<>) 2 , (28) 



where R = R/Hq, and i?o is its value today. Therefore from equation ([2 



fl„ = (29) 



and in general 



R = iT2 = — a ' y m ) 



R _2a- 2f3R + 7 - 3fi m e" 3A 

This relation is important as we know how R varies close to today. Taking the derivative of 
R with repect to time, plugging it into the 00 component of equation (|2"6"|) . we get the new 
Friedmann equation 

H 2 6n m e- 3N + f R R-f 



(31) 
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with H = H(N) is the Hubble factor H = a 1 da/dt, and 



= [2 e 3 ^/? 2 - 2a 7) -87 On] 2 
? 4[e 3Ar (/3 2 -2« 7 ) + 3 7 fiJ 2 



following 

The Friedmann equation fl3~Tl) together with its two iV-derivatives give us three equations 
for three unknowns. In contrast to the metric-based case, the Friedmann equation here does 
not involve any derivatives of the Hubble parameter, because R is independent of H. We 
end up with three non-linear equations involving Q m , 7i', Tl", to be solved simultaneously 
for a, j3 and 7. This time though, the equations are non-linear even in a, (3, and 7, so that 
there isn't a unique solution for these parameters. In general there is no telling how many 
real solutions there are. We solve the equations numerically in the range [-10,10] for the 
parameters. Listed in Appendix B are multiple solutions found at the mean values of the 
initial parameters. Next we consider one of these solutions. 

Linearization 

The equations of motion can be written as Fi(pj,q k ) = 0, with i = 1,2,3, where F\=0 
corresponds to the Friedmann equation, and -^2=0, and F 3 =0 to its two derivatives. Let 
Pj = (a,P, 7) and q k = (^l m , TC', TC"). The solutions listed in Appendix B are the pj which 
solve Fi(pj, q k ) = 0, where q k correspond to the mean values of the initial parameters. About 
any one solution we may linearize the equations of motion (for reasons explained in section 
II): 

(ij ~ Qj) ■ (33) 



dFi 

FiiPj, Qk) = w Fi(pj, q k ) + 



dp 



KPj ~Pj) + 



p.<i 



With Aij = dFi/dpj\p A and B ik = dFi/dq k \ M , one has 

Pj = -{A- 1 B) jk q k + pj + (A' 1 B) jk q k . (34) 

As an example, consider the solution [0.222102, 0.00488155, 4.12207e-05]. We can obtain 
the distributions of pj about this solution using the linearized equations together on the 
MCMC chains for the initial parameters. Results are shown in Fig 4, and for this solu- 
tion the matrix A~ x B is also given in Appendix B. Though the figure shows reasonably 
strong constraints on the modified gravity parameters, it should be remembered that in 
this formalism the data allow multiple such solutions. One may try to distinguish between 

11 



the solutions by imposing consistency checks derived from taking higher derivatives of the 
Friedmann equation. In that case we would use constraints on the Hubble parameter in 
more redshift bins, which would allow more freedom and hence more solutions overall to 
distinguish between (even though here we are only discussing solutions about one set of val- 
ues for the initial parameters). Hence this endevour will not be fruitful at this time. With 
much better data, the solutions may be tractable. 
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FIG. 4: Same as Fig 1, for f(R) gravity in the Palatini formalism. The constraints in the bottom 
panel are about one solution of the three non-linear equations of motion. There are other solutions, 
see text for further discussion. 

We also note in passing that when we constrain just two H n s, these turn out to be 
constrained to non-zero values (Fig 4). A negative Ti! implies that the universe is not super- 
accelerating today, and positive Ti" implies that Ti' is increasing (approaching zero) today. 
This is as expected. When more 7i"s are involved, as in the remaining sections, they of 
course become less constrained. 
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IV. f(R,R 2 GB ) METRIC-BASED THEORIES 



It would be interesting to look also at other generalizations which have been considered 



in the literature 



321 ]- 35]. These models have actions of the form 
1 



S 



16ttG 



d l x y /=g[R + f(R,R 2 GB )] + Sr, 



(35) 



where R GB = R 2 — 4i? Mi , R^ u + R fJiUpa R^ u P a is the Gauss-Bonnet combination. The equations 
of motion which follow can be written as 

(1 + F) G al3 - V a V p F + g aP OF - \g aP (f -FR-£R 2 GB ) 
-2RV a V P i + 2Rg aP Ui-A R a/3 □£ - 4 R {a ° T p) V CT V r £ 
- 4 g a(i R?° V P V^ + 8 R {a v V p) V„£ = 8nG T af3 , (36) 

where F = df/dR and £ = df/dR GB and the partial derivatives should be found treating 
R and R GB as independent variables. In a FRW background this becomes 



3 H 2 (1 + F + F') + \ (f - F R - £ R 2 GB ) + 12 ff 4 f = 87rG p . 



where 



p2 



24 - if 2 = 24 if 3 (if' + ff ) , 



(37) 



(38) 



and R has been already introduced in equations Qj. 

The Taylor-expansion of / about today up to second order in the scalars takes the form 



= a + (3(R-R ) + ^i(R-Ro) 2 

+ 72 (R - R ) (R gb - R GB0 ) + \ 73 (-R G b _ -^gbo) 



where 



'•■ = , /3 = F , 7! = if , 

#0 



<9F 
dR 



72 = ^df? 



^-"-GB 



(39) 



(40) 



with R = R/Hq and R GB = R gb /Hq. It should be noted that there is no linear term in 
the expansion for R GB , because such a term would give no contribution to the equations of 
motion, and the constant term in -RgbIo is considered absorbed in a. Then one has 



F = + Tl (R - Ro) + T2 (R 2 GB - R 2 GB0 ) 

H 2 t = l2 (R- Ro) + 73 (i^B - RIbo) ■ 



(41) 
(42) 
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Today equation (1371) becomes 

a = 6 [-1 + (1 + ft') (3 - 671 (4 ft' + ft' 2 + ft") 
-48 7 2 (4 ft' + 2H' 2 + ft") 

-96 73 (4 ft' + 3 ft' 2 + H") + Q m }. (43) 

Taking four derivatives of equation fl37j) allows us to solve for the five parameters in term of 
Q m and the the derivatives of ft (up to the 6th derivative). One derivative gives 

(3 = -1 + (2ft')" 1 [6(4ft' - 7ft' 2 - ft' 3 - 3ft" - 4ft'ft" - ft'") 7l ] 

+ (2ft')" 1 [48 (4ft' - 14ft' 2 - 6ft' 3 - 3ft" - 8ft'ft" - ft'")7 2 ] 

+ (2ft')- 1 [96(4ft' - 21ft' 2 - 15ft' 3 - 3ft" - 12ft'ft" - ft'") 73 ] - ^ (44) 

Further derivatives give 71, 72 and 73. The later two particularly involve extremely compli- 
cated expressions that cannot be written down here. 

When constraining 6 derivatives of the Hubble parameter from data, it is expected that 
the constraints will be poor. The top two panels of Fig 5 show the constraints on the initial 
parameters. For reasons discussed in section II, we will now linearize the relations for the 
modified gravity parameters about the mean likelihood values of the inital parameters. In 
this way we illustrate a space of possible solutions. (There would be a number of isolated 
solutions not represented in these results, but because these are isolated solutions they will 
also be suppressed in the full posterior. Better data should reduce the number of such 
solutions.) 

Linearization 

Linearizing the previous equations about the mean values of the the derivatives of ft and 
Q m one finds 

a = -4.52918 + 14.6151ft' - 3.28575ft" + -0.200232ft'" 

- 0.006872ft"" + 0.0004706ft'"" + 8.5558 x lO^ft""" + 18.3942ft m (45) 
(3 = -0.42128 + 4.91222ft'- 1.71918ft" + 0.1071ft'" 

- 0.003735ft"" + 0.000257ft'"" + 4.6535 x lO^ft""" + 5.39085fi m (46) 
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FIG. 5: Same as Fig 1, for f(R,RQ B ) gravity in the metric formalism. 



71 = -0.008276 - 0.05569ft' + 0.03140ft" - 0.002392ft'" 

+ 0.00004811ft"" - 6.2764 x 10" 6 ft""' - 1.03197 x 10~ 7 ft""" - 0.0705fi m (47) 

72 = 0.0006553 + 0.007342ft' - 0.005405ft" + 0.000428ft'" - 6.2403 x lO^ft"" 

+ 1.2073 x 10~ 6 ft""' + 1.9076 x 10~ 8 ft""" + 0.01237Q m (48) 
7s = 0.0003048 - 0.000411ft' + 0.0009694ft" - 0.00008026ft'" + 6.99103 x lO^ft"" 

- 2.5368 x lO^ft'"" - 3.8354 x 10~ 9 ft""" - 0.001939fi m (49) 

These relations together with the MCMC chains result in the constraints shown in the 
bottom panel of Fig 5. The parameters are all consistent with zero, the significance is in 
the order of the constraints. 

We have checked using the full equations that the 7's are in fact well constrained to be 
close to zero as indicated by the plot. This is because the relations contain combinations of 
the derivatives of ft in their denominators. Hence this is a consequence of the equations of 
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motion. This result is interesting. Again the H n s=0 is a singular point for the 7's (as well 
as for a and (5 the relations for which contain a division by TL'). Again there is a solution 
allowed by data that is a singular point in this theory, but overall (ie. for most allowed 
combinations of the derivatives of H) the parameters of the theory are well behaved, and 
the distributions shown encompass the main range of values for these parameters that are 
acceptable as per current data. 

V. DISCUSSION AND CONCLUSIONS 

f(R) and f(R, Rqb) models, by definition, introduce higher derivatives of the expansion 
rate into the equations of motion. This automatically requires data which can allow for 
precise measurements of derivatives of the Hubble parameter. This severe demand on ob- 
servations arises because currently we have no theory and no symmetry which can exactly 
predict the form of f{R). This problem is evidently similar to the task of finding the form for 
the potential of a quintessence-like field, or more generally to finding a dynamical alternative 
to a cosmological constant. On the other hand, we should not feel so complacent as to accept 
a tiny cosmological constant as the solution to the dark energy problem without giving due 
consideration to the daunting dynamical alternatives in the gravity sector. In order to get 
rid of the possibilities discussed in this paper of modifying gravity, one may introduce the 
axiom that only a linear combination of Lovelock terms can enter into the Lagrangian for 
gravity, although this axiom does not follow from any symmetry consideration. However, 
this would still not be enough to solve the cosmological constant problem, namely the zeroth 
Lovelock term, the cosmological constant that we would predict from QFT is far too large 
compared to the one needed by observations. f(R) models, or more generally f{R, Rqb) 
theories, thus remain a relevant consideration. 

We have used distance measurements from CMB, BAO and SNe la to place preliminary 
bounds on the parameters of f(R) models. f(R) was Taylor expanded about today, keeping 
terms upto second order in R. Equations of motion were found in both the metric and 
Palatini formalisms. These allow us to solve for the coefficients of the Taylor expansion of 
f{R). The equations of motion are non-linear, and we study the main space of possible 
solutions. The solutions are interesting in that they reflect the order of magnitude of the 
coefficients allowed by data. Under the metric formalism we find that df /dR\o is small 
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and negative over a substantial part of the allowed region, in which case in order for these 
models to have a consistent GR-like evolution at early times(see [23]), this quantity must 
have changed sign some time in the past. 

We have also analyzed a more general class of theories, f(R,R GB ). In order to modify 



gravity, without introducing spurious degrees of 
particles), one needs to use Lovelock scalars 
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reedom (such as extra ghost-like spin-2 



60| . In four dimensions the only such terms 
which give a non-zero contribution are a constant, R, and R GB - These models are more 
general and the constraints are weaker. However, as for the f(R) theories, we should make 
a serious attempt to shrink the allowed parameter space for such theories in different ways 
(using data and theoretical considerations). 

In addition to the analysis presented here, at the moment, other considerations (conver- 
gence to GR-like evolution at early times, instabilities) are still the most serious issues that 



these theories need to survive. For example, it is known that some models 
f(R, Rqb) do not have a FRW background consistent with GR at early times 



or f(R) and 
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6l|. Fur- 



thermore at high-redshift GR-like models may lead to unstable behaviour in the evolution of 



perturbations, because o 



speed of propagation 
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;he existence of either a tachionic mode or a mode with imaginary 
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621 ] . Such considerations will further constrain the parameter 



space of the models considered here. In this paper we have also not considered solar system 
type constraints, for reasons discussed in the Introduction. 

To conclude, we have used current data to place constraints on the first and second 
derivatives of general functions of the Ricci and the Gauss-Bonnet scalars, under different 
formalisms. These derivatives are important quantities for theoretical considerations relating 
to ghosts and instabilities. Next generation dark energy surveys 63j, |64], |65( can measure the 
cosmic expansion history much more precisely; this would dramatically shrink the presently 
allowed parameter space of the modified gravity models considered here. 
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APPENDIX A: INVERTING ff(z)'S TO THE DERIVATIVES OF H 



Consider the vector hi — Tii — 1 = H(zi) /Hq — 1. If the Hubble parameter is measured in 
n redshift bins, one can solve for n of its (present day) derivatives by assuming a truncated 
Taylor expansion for each of hi about today upto order n, which gives 



£ 

3=1 



dm 



5. 

2=0 J' 



(Al) 



This is a linear equation, hi = Aij , where = zj/j\, which can be inverted to find 
derivatives of H. (today). We then change variables from z to N, where N = — ln(l + z). 
One has 



dm 

~~dz} 



d 

d~N 



n 



so that 



dH 

dz 

dm 



dz 2 

dm 



dz 3 

dm 



dz 4 

dm 



dz 5 

dm 



dz 6 



-H'q 
nil! i nil 

H + H 

-(H Q 5) + 10 H ( o ] + 35 n { o ] + 50 H' Q ' + 24 H' ) 

nP + is + 85 + 225 + 274 nl + 120 n' 



In fact it can be shown that 



d k U 



dz k 



-l) k ^ akj Hp 

3=1 



(3) 



with 



afei = (k- 1)! 

akj = (k-1) Ofe-ij + Ofc-ij-i , 1 < j < k 

Q>kk = 1 • 



(A2) 

(A3) 
(A4) 
(A5) 
(A6) 
(A7) 
(A8) 

(A9) 



(A10) 
(All) 
(A12) 
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This then specifies the whole series of polynomials. 
Writing the linear relation as 

dm 



we have 



h i = A ij B jk nS e) = M^n^ , 



-1 nj(k) 



or 



= M kj hj 



For two bins case (zi = 0.7, and z 2 = 1.4 for current data) we get 



f -2.85714 0.714286^ 



M 



-1.22449 1.32653 



/ 



For three bins [z\ = 0.4667, z 2 = 0.9333, and z 3 = 1.4) then we get 

/ -6.4286 3.2143 -0.71429 \ 



M = 



-16.531 15.153 -3.8775 



\ 32.930 -22.369 3.2216 J 
for 4 bins equally spaced out to z — 1.4, we get 

( -11.429 8.5714 -3.8095 0.71428 ^ 



M 



-59.320 68.980 -34.286 6.7687 



-9.0962 55.802 -52.789 13.251 



509.11 -745.18 450.18 -91.607 



for 5 bins, we have 



/ 



M 



-17.857 17.857 -11.905 4.4647 -0.71429 



-145.83 209.61 -153.91 60.374 -9.9150 



-335.37 679.30 -630.53 276.88 -48.546 



1445.8 -2258.6 1643.0 -562.54 79.241 



2094.7 -6287.7 7809.3 -4285.9 838.55 



\ 
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and, finally, for 6 bins, we have 



/ 


-19.3 


16.109 


-7.193 


0.03758 


1.271 


-0.3546 \ 




-179.35 


218.78 


-108.04 


0.67106 


20.4867 


-5.823 




-500.66 


833.8 


-523.69 


4.9750 


118.16 


-35.142 




1837.10 


-2403.5 


1168.7 


4.745 


-149.15 


32.4457 




4187.9 


-9101.2 


7205.2 


-112.35 


-2310.5 


732.79 


V 


-54953 


91546 


-59415 


191.52 


14897 


-4190.4 / 



(A20) 



With the present day derivatives of the Hubble parameter thus derived from model indepen- 
dent measurements of the Hubble parameter in linear redshift bins from cosmological data, 
we can then obtain constraints on parameters of modified gravity models, using appropriate 
equations of motion. This is done in sections II, III and IV of this paper. 



APPENDIX B: APPENDIX B: SOLUTIONS FOR THE PALATINI FORMALISM 



We could find the following distinct solutions of the equations Fi(pj,qk) = for 
a, 0, and 7, about the mean values for q k = (a m .H',H"): [-0.857943, -0.428971, - 
0.14299], [-1.46637, -0.422284, -0.140728], [0.399627, 0.824305, 0.852729], [3.78266, 0.748323, 
-0.00792803], [0.222102, 0.00488155, 4.12207e-05], [-0.739164, -0.00165984, 0.0000280094], 
[-0.866737, -0.00527848, 0.000315044], [0.252698, -0.00306054, 0.0000152404], [-0.29676, 
0.00333869, 0.0000390608], [-0.811148, -0.00218394, 0.0000614933], [-0.410838, -0.0124848, 
-0.0000849403], [-0.132685, -0.00276524, 0.0000205499], and [0.00810481, -0.00504213, - 
0.00002751]. 

For the solution [0.222102, 0.00488155, 4.12207e-05], the A~ l B of equation (El) is 



/ 



A- l B 



8.7943 -0.0077692 0.000655508 
11.162 -0.015052 0.0012888 
0.18687 -0.00025239 0.000021581 



\ 



:bi) 
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